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Simulated aeroelastic responses of a nonlinear pitch and plunge apparatus are analyzed 
using various statistical signal processing techniques including higher-order spectral meth- 
ods. A MATLAB version of the Nonlinear Aeroelastic Testbed Apparatus (NATA) at the 
Texas A&M University is used to generate various aeroelastic response data including limit 
cycle oscillations (LCO). Traditional and higher-order spectral (HOS) methods are applied 
to the simulated aeroelastic responses. Higher-order spectral methods are used to identify 
critical signatures that indicate the transition from linear to nonlinear (LCO) aeroelastic 
behavior. 


Introduction 

The study of nonlinear aeroelastic phenomena, in particular that of limit cycle oscillations (LCO), has 
become an important field of study in recent years. 1 A major motivating factor for this increased interest is 
the fact that several operational aircraft continue to experience LCO within their flight envelope, often leading 
to a flight restriction or performance degradation. 2,3 Several researchers have been using the Nonlinear 
Aeroelastic Testbed Apparatus (NATA) at Texas A&M University’s Low-Speed Wind Tunnel to investigate 
fundamental aspects of LCO behavior. 4-9 

Recently, higher-order spectral methods have been applied to experimental flutter measurements from 
wind-tunnel data 10 and flight test data. 11 Analyses using higher-order spectral methods are proving beneficial 
to improved understanding of nonlinear aeroelastic phenomena. Higher-order spectral methods enable the 
visualization of the transfer of energy from one frequency to another, a hallmark of nonlinear phenomena 
that is not visible using traditional spectral methods. 

The paper is organized as follows. A detailed description of the NATA and its MATLAB model are 
provided, including the equations of motion that define the dynamics of the system. A description of higher- 
order spectral (HOS) methods is provided, including an illustrative example. Simulated data generated 
using a variety of control surface inputs are then analyzed using traditional and HOS methods. The issue 
of the Gaussianess of an input/output combination for identifying regions of linear and nonlinear behavior 
is presented and applied. 

Nonlinear Aeroelastic Testbed Apparatus (NATA) 

The NATA in the Texas A&M University’s 2 ft. x 3 ft. Low Speed Wind Tunnel is used to generate 
free and controlled, linear and nonlinear aeroelastic transients at various conditions. Based on a simple 
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two-degree-of-freedom concept with structurally nonlinear components, a rigid wing mounted on the NATA 
can exhibit nonlinear aeroelastic behavior, including LCO phenomena. The pitch and plunge motions of the 
NATA are enabled by springs that can be tuned to a wide variety of nonlinear stiffness characteristics. The 
NATA is shown in figure 1. 



Figure 1. Schematic of the Nonlinear Aeroelastic Testbed Apparatus (NATA). 

A rigid wing attached to the NATA, as it would be mounted in the Low Speed Wind Tunnel, is shown 
in figure 2. The nonlinear structural response is governed by a pair of cams that can be tuned to model 
various types of polynomial stiffness functions. Additional details regarding the NATA are available in the 
references. 4-8 


Analytical Model of NATA 


A MATLAB model of the NATA/rigid wing system was generated and is briefly discussed in this section. 
The equations of motion for the system are presented as Eq. 1 
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m w mass of the wing 
x a static unbalance 


2 of 20 


American Institute of Aeronautics and Astronautics 



Plunge 



Figure 2. Schematic of the Nonlinear Aeroelastic Testbed Apparatus (NATA) with rigid airfoil and control 
surface. 

b semichord 

I a mass moment of inertia about the elastic axis 

h plunge degree of freedom 

a pitch degree of freedom 
Ch plunge damping coefficient 
c a pitch damping coefficient 
kh plunge stiffness coefficient 
k a pitch stiffness coefficient 
L lift 

M moment 


The lift (L) and moment (M) on the right-hand side of Eq. 1 are defined using quasi-steady aerodynamics 
as presented in Eq. 2 and Eq. 3. 

L = pU 2 bcia ^ + (d - oj b^j + pU 2 bci 0 (3 (2) 

M = pU 2 b 2 c ma ^ + Q - oj b^j + pU 2 b 2 c m0 (3 (3) 

where 

p density 

U freestream velocity 

ci a lift due to angle of attack 

cip lift due to control surface deflection 

/? control surface deflection 

c ma moment due to angle of attack 

c m/ 3 moment due to control surface deflection 
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Higher-Order Spectra 


The primary benefit of higher-order spectra (HOS), also known as higher-order frequency response func- 
tions, is that they provide information regarding the interaction of frequencies due to a nonlinear process. 
For example, bispectra have been used in the study of grid-generated turbulence to identify the nonlinear ex- 
change of energy from one frequency to another (related to the turbulent cascade phenomenon). Traditional, 
or standard, concepts, by definition, cannot provide this type of information. In addition, higher-order spec- 
tra are the frequency-domain version of the Volterra series. For details regarding this relationship, the reader 
is referred to the recent article by Silva . 12 Some very interesting and fundamental applications using the 
frequency-domain Volterra theory 13, 14 and experimental applications of Volterra methods 15, 16 are providing 
new “windows” on the world of nonlinear aeroelasticity. 

In the recent work by Hajj and Silva , 17, 18 the aerodynamic and structural aspects of the flutter phe- 
nomenon of a wind-tunnel model are determined via a frequency domain analysis based on a hierarchy of 
spectral moments. The power spectrum is used to determine the distribution of power among the frequency 
components in the pressure, strain and acceleration data. The cross-power spectrum, linear coherence, and 
phase relation of the same frequency components between different signals are used to characterize the 
bending and torsion characteristics of the model. The nonlinear aspects of the aerodynamic loading are 
determined from estimates of higher-order spectral moments, namely, the auto- and cross-bispectrum. 

For a discrete, stationary, real-valued, zero-mean process, the auto-bispectrum is estimated as 19 

M 

B xxx [h ,h] = ^Yl X T ] [h + 1 2 }K {k) [h]X* T (k) [h] (4) 

where Xj?\l] is the Discrete Fourier Transform of the k th ensemble of the time series x(t) taken over a 
time, T, and M is the number of these ensembles. The auto-bispectrum of a signal is a two-dimensional 
function of frequency and is generally complex- valued. In averaging over many ensembles, the magnitude of 
the auto-bispectrum will be determined by the presence (or absence) of a phase relationship among sets of 
the frequency components at Zi , Z2, and l\ + I2. If there is a random phase relationship among these three 
components, the auto-bispectrum will average to a very small value. Should a phase relationship exist among 
these frequency components, the corresponding auto-bispectrum will have a large magnitude . 20 Because a 
quadratic nonlinear interaction between two frequency components, li and I2, yields a phase relation between 
them and their summed component, Zi H- Z2? the auto-bispectrum can be used to detect a quadratic coupling 
or interaction among different frequency components of a signal. The level of such coupling in a signal 
can then be associated with a normalized quantity of the auto-bispectrum, called the auto-bicoherence, and 
defined as 
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By the Schwarz inequality, the value of b^ xx [lu 1 2] varies between zero and one. If no phase relationship 
exists among the frequency components at Zi , Z2, and l\ + Z2, the value of the auto-bicoherence will be at or 
near zero (due to averaging effects). If a phase relationship does exist among the frequency components at 
Zi , I2 1 and li + /2> then the value of the auto-bicoherence will be near unity. Values of the auto-bicoherence 
between zero and one indicate partial quadratic coupling. 

For systems where multiple signals are considered, detection of nonlinearities can be achieved by using 
the cross-spectral moments. For two signals x(t) and y(t), their cross-bispectrum is estimated as 


B yxx [hM = T E W + h]X* T {k) [h]X* T (k) [h] (6) 

1 k = 1 

where Xj?\l\ and Y^\l\ are the Discrete Fourier Transforms of the k th ensemble of the time series x(t) 
and y(t), respectively, over a time, T. The cross-bispectrum provides a measure of the nonlinear relationship 
among the frequency components at li and I2 in x(t) and their summed frequency component, li + Z2, in 
y(t). Similar to the auto-bispectrum, the cross-bispectrum of signals x(t) and y(t) is a two-dimensional 
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function in frequency and is generally complex- valued. In averaging over many ensembles, the magnitude 
of the cross-bispectrum will also be determined by the presence, or absence, of a phase relationship among 
sets of the frequency components at Zi , Z 2 , and l\ + h- If there is a random phase relationship among the 
three components, the cross-bispectrum will average to a very small value. Should a phase relationship exist 
among these frequency components, the corresponding cross-bispectral value will have a large magnitude. 
The cross-bispectrum is thus able to detect nonlinear phase coupling among different frequency components 
in two signals because of its phase-preserving effect. 

Similar to defining the auto-bicoherence, one can define a normalized cross-bispectrum to quantify the 
level of quadratic coupling in two signals. This normalized value is called the cross-bicoherence and is defined 
as 
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If no phase relationship exists among the frequency components at Zi , I 2 in x(t) and the frequency 
component at 1% + Z 2 in 2/(t), the value of the cross-bicoherence will be at or near zero. If a phase relationship 
does exist among these frequency components, the value of the cross-bicoherence will be near unity. Values 
of cross-bicoherence between zero and one indicate partial quadratic coupling. A digital procedure for 
computing the auto and cross-bicoherence is given by Kim and Powers 19 and is summarized by Hajj et al. 21 

An important property of HOS functions involves the Gaussianess of a signal. The auto-bicoherence of a 
Gaussian signal is identically zero. This is an important theoretical result that is described in the literature. 
In addition, it is well-known that if a Gaussian signal is input to a linear system, then the output will be 
Gaussian as well. Therefore, computation of the auto-bicoherence for a Gaussian output will yield a value at 
or near zero, indicative of a linear process. If a Gaussian signal is applied to a nonlinear system, the output 
will be non-Gaussian. As a result, the auto-bicoherence for the output of the nonlinear system will yield 
a non-zero value. This characteristic of HOS methods is exploited for the present application to identify 
regions where the behavior is transitioning from linear to nonlinear dynamics. 

Using a linear (Eq. 8) and a nonlinear differential equation (van der Pol oscillator, Eq. 9), the effect of 
nonlinear processes on input Gaussian distributions can be demonstrated. A Gaussian signal is input to 
the linear and nonlinear differential equations. The histogram of the response for each system (linear and 
nonlinear) to the Gaussian input signal is presented in the following figures. 

Presented in Figure 3 are the input and output histograms for the linear system. As can be seen, the 
linear system preserves the Gaussianess of the input signal. Presented in Figure 4 are the input and output 
histograms for the nonlinear system. For the nonlinear case, the Gaussianess of the input signal has not 
been preserved and has, instead, been transformed into a non-Gaussian distribution radically different from 
the input distribution. This is a well-known characteristic of nonlinear systems. It is worthwhile to mention 
that the response of the van der Pol oscillator (the nonlinear ODE) is a limit cycle oscillation (LCO). It 
appears that the distribution for other LCO phenomena may have a similar structure. 

HOS methods are very sensitive to variations in Gaussianess and, as will be shown, are valuable in 
identifying deviations from Gaussian behavior. Deviations from Gaussian behavior, as shown above, indicate 
the existence of a nonlinear process. 


Vi = 2 / 2 , 2/2 = ~V2 ~yi+ input 


(8) 


m = U2, m = - m (! - 2 / i ) 2 2/2 - m + input 


(9) 


Example 2 

A second example is now presented in order to explain some of the concepts associated with HOS. The 
example 22 consists of a linear (y) and a nonlinear (z) time series, both containing Gaussian noise (d). The 
equations for these time series are 

x(k) = sin(27r(/i)fc + |) + sin(27r(/ 2 )fc + |) (10) 
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Input to Linear Differential Equation 



Output from Linear Differential Equation 


Figure 3. Histograms of the input and output for the linear system due to a Gaussian input. 



Input to Nonlinear Differential Equation 



Output from Nonlinear Differential Equation 


Figure 4. Histograms of the input and output for the nonlinear system due to a Gaussian input. 


6 of 20 


American Institute of Aeronautics and Astronautics 


y(k) = x(k) + d(k) 


( 11 ) 


z{k) = x(k) + 0.05x 2 (&) + d(k) (12) 

where /i = 0.12 Hz and = 0.30 Hz. 

Figures 5 and 6 each contain the time series and the associated magnitude of the Fourier transforms for 
the linear and nonlinear time series, respectively. The dominant frequencies (in this case 0.12 and 0.30 Hz) 
are clearly visible in both figures. The frequency content for the nonlinear time series indicates the existence 
of additional frequencies. The nature of these frequencies, whether or not these frequencies are random or 
the result of a nonlinear coupling process, cannot be discerned from this analysis. 




Figure 5. Linear time series and magnitude of frequency response. 


The power spectral densities (PSD) for the linear and nonlinear time series are presented as Figure 7. 
Here again, the existence of additional frequencies in the PSD of the nonlinear time series (as compared to 
the linear time series) is obvious. However, the PSD information cannot be used to discern if these additional 
frequencies are random or the result of a nonlinear coupling process. 

In order to determine if these additional frequencies are random or the result of a nonlinear coupling 
process, the HOS for these time series must be computed. The magnitude of the auto-bicoherence for the 
linear time series, y, is presented as Figure 8. The two frequency axes correspond to the two frequency 
indices in Eq. 5. If a significant peak is observed, this implies that the sum of those two frequencies is the 
result of a quadratic (nonlinear) coupling. It can be seen that the magnitude of the auto-bicoherence for 
the linear time series is quite low in value (compared to unity) and that there exist no dominant peaks that 
would be indicative of a nonlinear coupling. 

The magnitude of the auto-bicoherence for the nonlinear time series is presented as Figure 9. In contrast 
to the auto-bicoherence for the linear time series, the auto-bicoherence for the nonlinear time series has 
several peaks at or close to unity, indicative of nonlinear (quadratic) coupling at the x- and y-axis frequencies 
indicated. 

Typically, to enhance visualization and interpretation of the auto- or cross-bicoherence , contour plots are 
viewed. The contour plot for Figure 9 is presented as Figure 10. The contour plot presents the frequencies 
that have coupled quadratically to generate the new frequencies that were visible in the PSD. In addition, 
the symmetry associated with the computation of the auto-bicoherence, in this case, is also evident in the 
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Time Steps 



Figure 6. Nonlinear time series and magnitude of frequency response. 



Figure 7. Power spectral density (PSD) functions for linear and nonlinear time series. 
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Figure 9. Auto-bicoherence for the nonlinear time series. 
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contour plot. The bicoherence basically computes the correlation between two frequencies and the sum of 
those two frequencies. For example, the magnitude of the auto-bicoherence when the first frequency (x-axis) 
is 2 Hz and the second frequency (y-axis) is 5 Hz will be the same as the magnitude of the auto-bicoherence 
when the first frequency is 5 Hz and the second frequency is 2 Hz since their sum is the same. 


i r- 
0.8 - 
0.6 - 
0.4 - 


n 0-2 - 



- 0.4 - ■ .... 

- 0.6 - 
- 0.8 - 

l l l l l l l l l 

-1 - 0.8 - 0.6 - 0.4 - 0.2 0 0.2 0.4 0.6 0.8 1 

Frequency, Hz 

Figure 10. Contour plot of the auto-bicoherence for the nonlinear time series. 


This symmetry is presented in Figure 11. As Figure 11 indicates, knowledge of Region I (two such regions 
in the first quadrant) and Region II (two such regions in the fourth quadrant) is sufficient to completely 
define the remainder of the quadrants due to symmetry considerations. Regions I and II indicate that the 
original frequencies have been quadratically-coupled, resulting in new frequencies that are the sum (positive 
and negative) of the original frequencies. The new frequencies include 0.18, 0.24, and 0.42. The limitation 
of Region I to the triangle shown is due to the fact that the summation of two frequencies cannot exceed 
the Nyquist frequency (0.5, in this case). Therefore, the combination of the second frequency (0.30 Hz) with 
itself is not included since 0.60 Hz would be greater than the Nyquist frequency. 

For subsequent results, although the auto-bicoherence is computed for all frequencies (positive and neg- 
ative), Region I in the first quadrant will be of primary interest for the sake of simplicity. This will be 
sufficient to demonstrate the applicability and value of HOS methods to aeroelastic data. 

Results 

In this section, aeroelastic transients that range from a stable response to an LCO response, generated 
in MATLAB using the analytical model of the NATA, are presented. Then, in order to investigate the level 
of nonlinearity for each transient, Gaussian inputs will be applied to the system via the control surface. The 
output responses will be analyzed for their level of Gaussianity. Deviations from a Gaussian distribution can 
be interpreted as deviations from linear behavior. The auto-bicoherence is also presented to demonstrate the 
ability of HOS methods to discern Gaussianess from non-Gaussianess. For simplicity, the results presented 
are limited to the pitch motion. 

Figure 12 presents the aeroelastic transient in pitch for a condition where a stable aeroelastic response 
is generated at a small velocity (U=5 m/sec) with small initial conditions as perturbations and no control 
surface input. The aeroelastic transient is clearly stable and the frequencies of the system are evident in the 
plot of the magnitude of the frequency response function of the pitch motion. 

Figure 13 presents the aeroelastic transient in pitch for a condition where an LCO has been reached at a 

10 of 20 


American Institute of Aeronautics and Astronautics 



Frequency, Hz 


Figure 11. 


Contour plot of the auto-bicoherence for the nonlinear time series with symmetry regions identified. 



Frequency, Hz 


Figure 12. 


Stable aeroelastic response in pitch at U=5 m/sec: Time domain and frequency domain. 
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velocity of U=20 m/sec. The frequency content is exhibiting greater complexity than the previous aeroelastic 
transient. 

To determine the level of nonlinearity that may be present starting from the stable transient (U=5 m/sec) 
all the way to the LCO transient (U=20 m/sec), a Gaussian input is applied via the control surface and the 
resulting nature of the output transient (Gaussian or non- Gaussian) is determined. This approach serves to 
demonstrate the sensitivity of HOS methods to non- Gaussian signals. 



10 10 

Frequency, Hz 


Figure 13. Limit cycle oscillation (LCO) aeroelastic response in pitch at U=20 m/sec: Time and frequency 
domain. 


A Gaussian input is generated and applied to the aeroelastic system via the control surface. The histogram 
of the input is presented as Figure 14, where the Gaussian nature of the signal is evident. The auto- 
bicoherence for the Gaussian input is presented as Figure 15. The small values of the auto-bicoherence for 
the Gaussian input are as expected for a signal with a Gaussian distribution. 

The resultant aeroelastic response (at U=5 m/sec) and the corresponding magnitude of the frequency 
response are presented in Figure 16. The effect of the Gaussian input is evident in the transient as well as 
the frequency content when Figure 16 is compared with Figure 12. However, neither of these two functions 
can be used to indicate the presence or level of any nonlinear interactions. The histogram of the output 
response (pitch, in this case) is presented as Figure . 

The histogram for the output response indicates a noticeable deviation from the theoretical Gaussian 
distribution , indicating that the Gaussian input has been processed through some level of nonlinear inter- 
actions. 

The auto-bicoherence for this aeroelastic response, presented as Figure 18, confirms the existence of 
nonlinear coupling at frequencies close to the LCO frequency (3 Hz). 

The velocity is increased to U=10 m/sec and the same Gaussian input is applied. The aeroelastic 
response, along with the magnitude of the frequency response, is presented in Figure 19. 

Once again, the information presented in Figure 19 cannot be used to discern if any nonlinear interactions 
are taking place at this condition. However, the histogram of the output, presented as Figure 20, indicates 
an increased level of deviation from Gaussianess. 

The auto-bicoherence for this condition is presented as Figure 21, with increasing levels of nonlinear 
interaction near the frequency of the LCO. The important point to be made is that both, the histograms 
and the auto-bicoherence, are indicating the presence of nonlinear phenomena prior to the visible onset of 
the LCO at U=20 m/sec. 
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Control Surface Deflection, Degrees 


Figure 14. 


Histogram of the Gaussian input and theoretical Gaussian distribution (red). 
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Figure 15. Auto-bicoherence of the Gaussian control surface input. 
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Frequency, Hz 


Figure 16. Stable aeroelastic response in pitch at U=5 m/sec with Gaussian control surface input: Time and 
frequency domain. 
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Figure 18. Auto-bicoherence of stable aeroelastic response in pitch at U=5 m/sec due to Gaussian control 
surface input. 




Figure 19. Stable aeroelastic response in pitch at U=10 m/sec with Gaussian control surface input: Time and 
frequency domain. 
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Figure 20. Histogram of stable aeroelastic response in pitch at U=10 m/sec due to Gaussian control surface 
input. 


Auto-Bicoherence 



Figure 21. Auto-bicoherence of stable aeroelastic response in pitch at U=10 m/sec due to Gaussian control 
surface input. 
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The same Gaussian control surface input was applied for the U=20 m/sec condition at which LCO is 
encountered. The resulting LCO aeroelastic transient and the corresponding magnitude of the frequency 
response in pitch for that conditions is presented as Figure 22. The histogram of the pitch response is 
presented as Figure 23, where the deviation from Gaussianess is obvious. The resulting histogram of the 
LCO pitch response is clearly non-Gaussian with a bi-modal character. Figure 24 presents the corresponding 
aut-bicoherence for this condition. The auto-bicoherence for this LCO condition captures the nonlinear 
interactions that occur in the vicinity of the LCO frequency with greater resolution than the previous 
conditions. 


0.2 



Time, sec 



Figure 22. LCO aeroelastic response in pitch at U=20 m/sec with Gaussian control surface input: Time and 
frequency domain. 

Therefore, using this type of information (input/output histograms and auto-bicoherence) may enable 
the identification of precipitous nonlinear behavior prior to the actual onset during wind-tunnel or flight 
testing. This capability would provide greater insight into the testing process as well as enhanced insight 
regarding the nature of the nonlinear aeroelastic response. 

Conclusion 

Simulated responses from a MATLAB model of the Texas A&M University’s Nonlinear Aeroelastic 
Testbed Apparatus (NATA) with a rigid wing were presented. Aeroelastic transients generated included 
stable aeroelastic responses as well as nonlinear responses, or limit cycle oscillations (LCO). A Gaussian 
input was applied to the system via the control surface in order to determine the deviation from Gaussian 
behavior (i.e., deviation from linear behavior) of the aeroelastic system. Higher-order spectra were computed 
and served to indicate the existence of nonlinear interactions prior to the LCO event. 
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600 



Fitch, degrees 


Figure 23. Histogram of stable aeroelastic response in pitch at U=10 m/sec due to Gaussian control surface 
input. 



Figure 24. Auto-bicoherence of stable aeroelastic response in pitch at U=10 m/sec due to Gaussian control 
surface input. 
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